function out = func_LL(theta,choices,choice_set,covariates)

% Likelihood function 
%
% Inputs:
% - type: 
%   'LL': log-likelihood 
%   'GRADc': vector of individual likelihoods when extra term is included (needed for computation of OPG)
% - theta: parameter estimate at which to evaluate the function
% - choice_set: choice set matrix (1 if school in choice set, 0 otherwise)
% - choices: matrix of choices (1 if school is chosen, 0 otherwise)
%
%  Output:
% - out: 
%   log-likelihood function evaluated at theta

% ----------------------- %
% 1/ CHOICE PROBABILITIES %
% ----------------------- %

% # of parameters
% nb_params = size(theta,1);

% # of (exploded) observations
Nobs = size(choices,1);

% # of schools
J = size(choice_set,2);

% theta(1:end-1) = theta(1:end-1)./exp(theta(end));
% coeff_dis_est = -1./exp(theta(end));
    
%disp(mat2str(theta,3));

% V_ij: Fixed component of utility (N x J matrix)
% V_ij = repmat([0, theta(1:J-1)'], [Nobs 1]);
V_ij = repmat(theta(1)*[1:J], [Nobs 1]);

replication_factor = size(choices,1)/size(covariates,1);

% if nb_params>J-1
%     for pp=J:nb_params          
%         V_ij = V_ij + theta(pp).*repmat(covariates(:,:,pp+1),replication_factor,1);
%     end       
% end

V_ij = V_ij + theta(2).*repmat(covariates(:,:,J+1),replication_factor,1) ...
    + theta(3)*repmat(covariates(:,:,J+2),replication_factor,1) ...
    + theta(4)*repmat(covariates(:,:,J+3),replication_factor,1);

% % If model includes distance to school
% if isempty(distances)==0
%     % Distance matrix has to be duplicated to match the (exploded) choice set: 
%     replication_factor = Nobs/IA;
%     distance_matrix = repmat(distances, [replication_factor 1]); 
%     % Fixed component of utility evaluated at theta 
%     V_ij = V_ij + theta(J).*distance_matrix;
% end
% 
% % If model includes own score x school score
% if isempty(sch_stu_scores)==0
%     % Score matrix has to be duplicated to match the (exploded) choice set: 
%     replication_factor2 = Nobs/IA;
%     scores_matrix = repmat(sch_stu_scores, [replication_factor2 1]); 
%     % Fixed component of utility evaluated at theta   
%     V_ij = V_ij + theta(J+1).*scores_matrix;
% end

% exp(V_ij)
expV = exp(V_ij);

% If the model includes only school f.e., use version below (faster)
% expV = repmat(exp([0, theta']), [Nobs 1]);

% Probability that each school is preferred in the choice set (IA x J matrix)
prob = expV.*choice_set./repmat((expV.*choice_set)*ones(J,1), [1 J]);

% Log of choice probability for chosen school
log_prob = log(prob(choices==1));

% prob_choice = prob(choices==1);
% verif = (choices==1);
% prob_verif = sum((prob.*verif)',1);
% test = prob_verif';


% --------- %
% 3/ OUTPUT %
% --------- %

% Output: log-likelihood
out = -sum(log_prob);

   